# QJPS Replication
# Figures 1 and 2
# Updated by David K. Park

# Set Working Directory HERE
# setwd("C:/Research/Blue Red States/QJPS/Data/Replication/Replication1-2")
setwd("C:/Andy/Dropbox/ReplicationProject/Gelman") #desktop
setwd("C:/Duke/Dropbox/ReplicationProject/Gelman") #laptop
setwd("C:/Users/aob5/Dropbox/ReplicationProject/Gelman") #Bunche
# load foreign package to read in STATA file
library(foreign)

fips.cbs <- read.dta("fips.icpsr.cbs.naes.dta")

# Figure 1 Regressions
# state-level
st <- read.table("st.dat")

st.all <- merge(st, fips.cbs, by.x="stateabb", by.y="st")
st.s <- subset(st.all, st.all$regioncbs==3)
st.ns <- subset(st.all, st.all$regioncbs!=3)

p.st.yr <- seq(min(st$year), max(st$year),4)
n.st.yr <- length(p.st.yr)

col.st.names <- c("st.int", "z.st.inc.pop")
  n.st.coef <- length(col.st.names)

B.st.lm <- array(NA, c(n.st.yr, n.st.coef))
  colnames(B.st.lm) <- col.st.names
SE.st.lm <- array(NA, c(n.st.yr, n.st.coef))
  colnames(SE.st.lm) <- col.st.names

B.st.s.lm <- array(NA, c(n.st.yr, n.st.coef))
  colnames(B.st.s.lm) <- col.st.names
SE.st.s.lm <- array(NA, c(n.st.yr, n.st.coef))
  colnames(SE.st.s.lm) <- col.st.names

B.st.ns.lm <- array(NA, c(n.st.yr, n.st.coef))
  colnames(B.st.ns.lm) <- col.st.names
SE.st.ns.lm <- array(NA, c(n.st.yr, n.st.coef))
  colnames(SE.st.ns.lm) <- col.st.names

for (i in 1:n.st.yr){
    st.temp <- subset(st, st$year==p.st.yr[i])
    st.reg <- summary(lm(st.temp$st.repshare ~ st.temp$z.st.inc.pop))
    B.st.lm[i,] <- st.reg$coef[,1]
    SE.st.lm[i,] <- st.reg$coef[,2]
}

for (i in 1:n.st.yr){
    st.s.temp <- subset(st.s, st.s$year==p.st.yr[i])
    st.s.reg <- summary(lm(st.s.temp$st.repshare ~ st.s.temp$z.st.inc.pop))
    B.st.s.lm[i,] <- st.s.reg$coef[,1]
    SE.st.s.lm[i,] <- st.s.reg$coef[,2]
}

for (i in 1:n.st.yr){
    st.ns.temp <- subset(st.ns, st.ns$year==p.st.yr[i])
    st.ns.reg <- summary(lm(st.ns.temp$st.repshare ~ st.ns.temp$z.st.inc.pop))
    B.st.ns.lm[i,] <- st.ns.reg$coef[,1]
    SE.st.ns.lm[i,] <- st.ns.reg$coef[,2]
}

B.st.lm <- data.frame(B.st.lm)
SE.st.lm <- data.frame(SE.st.lm)

B.st.s.lm <- data.frame(B.st.s.lm)
SE.st.s.lm <- data.frame(SE.st.s.lm)

B.st.ns.lm <- data.frame(B.st.ns.lm)
SE.st.ns.lm <- data.frame(SE.st.ns.lm)

#  Figure 1 Plot
windows(width=8, height=2.5)
par(mfrow=c(1,3))
n.st <- NROW(B.st.lm)

plot(0, 0, type='n', ylim=c(-.3, .3), xlim=c(min(st$year), max(st$year)),
  cex.lab=1.15, xaxt="n", yaxt="n",  main="All States", cex.main=1, xlab="Year", ylab="Income Coefficient")
  axis(side=1, at=c(1960, 1980, 2000), cex.axis=1.1)
  axis(side=2, at=c(-.2, 0, .2), cex.axis=1.1)
  abline(h=0, col="gray")
    for (i in 1:n.st) {
    points(p.st.yr[i], B.st.lm[i,2], type='p', pch=19)
    segments(p.st.yr[i], (B.st.lm[i,2]-SE.st.lm[i,2]), p.st.yr[i], (B.st.lm[i,2]+SE.st.lm[i,2]))
  }

plot(0, 0, type='n', ylim=c(-.3, .3), xlim=c(min(st$year), max(st$year)),
  cex.lab=1.15, xaxt="n", yaxt="n",  main="Southern States", cex.main=1, xlab="Year", ylab="Income Coefficient")
  axis(side=1, at=c(1960, 1980, 2000), cex.axis=1.1)
  axis(side=2, at=c(-.2, 0, .2), cex.axis=1.1)
  abline(h=0, col="gray")
    for (i in 1:n.st) {
    points(p.st.yr[i], B.st.s.lm[i,2], type='p', pch=19)
    segments(p.st.yr[i], (B.st.s.lm[i,2]-SE.st.s.lm[i,2]), p.st.yr[i], (B.st.s.lm[i,2]+SE.st.s.lm[i,2]))
  }

plot(0, 0, type='n', ylim=c(-.3, .3), xlim=c(min(st$year), max(st$year)),
  cex.lab=1.15, xaxt="n", yaxt="n",  main="Non-Southern States", cex.main=1, xlab="Year", ylab="Income Coefficient")
  axis(side=1, at=c(1960, 1980, 2000), cex.axis=1.1)
  axis(side=2, at=c(-.2, 0, .2), cex.axis=1.1)
  abline(h=0, col="gray")
    for (i in 1:n.st) {
    points(p.st.yr[i], B.st.ns.lm[i,2], type='p', pch=19)
    segments(p.st.yr[i], (B.st.ns.lm[i,2]-SE.st.ns.lm[i,2]), p.st.yr[i], (B.st.ns.lm[i,2]+SE.st.ns.lm[i,2]))
  }

# Figure 1
savePlot(filename="st.yr", type=c("pdf"), device=dev.cur())
savePlot(filename="st.yr", type=c("ps"), device=dev.cur())


# individual level
anes <- read.table("recode.anes.qjps.txt")
anes.s <- subset(anes, anes$regioncbs==3)
anes.ns <- subset(anes, anes$regioncbs!=3)

# table(anes$relatt, anes$yr)
# table(anes$relatt2, anes$yr)
p.yr <- seq(min(anes$yr), max(anes$yr),4)
n.yr <- length(p.yr)
col.names <- c("int", "z.inc")
n.coef <- length(col.names)

B.glm <- array(NA, c(n.yr, n.coef))
  colnames(B.glm) <- col.names
SE.glm <- array(NA, c(n.yr, n.coef))
  colnames(SE.glm) <- col.names

B.s.glm <- array(NA, c(n.yr, n.coef))
  colnames(B.s.glm) <- col.names
SE.s.glm <- array(NA, c(n.yr, n.coef))
  colnames(SE.s.glm) <- col.names

B.ns.glm <- array(NA, c(n.yr, n.coef))
  colnames(B.ns.glm) <- col.names
SE.ns.glm <- array(NA, c(n.yr, n.coef))
  colnames(SE.ns.glm) <- col.names

for (i in 1:n.yr){
    temp <- subset(anes, anes$yr==p.yr[i])
    reg <- summary(glm(temp$vote ~ temp$z.inc, family=binomial(link=logit)))
    B.glm[i,] <- reg$coef[,1]
    SE.glm[i,] <- reg$coef[,2]
}

for (i in 1:n.yr){
    temp.s <- subset(anes.s, anes.s$yr==p.yr[i])
    reg.s <- summary(glm(temp.s$vote ~ temp.s$z.inc, family=binomial(link=logit)))
    B.s.glm[i,] <- reg.s$coef[,1]
    SE.s.glm[i,] <- reg.s$coef[,2]
}

for (i in 1:n.yr){
    temp.ns <- subset(anes.ns, anes.ns$yr==p.yr[i])
    reg.ns <- summary(glm(temp.ns$vote ~ temp.ns$z.inc, family=binomial(link=logit)))
    B.ns.glm[i,] <- reg.ns$coef[,1]
    SE.ns.glm[i,] <- reg.ns$coef[,2]
}

B.glm <- data.frame(B.glm)
SE.glm <- data.frame(SE.glm)

B.s.glm <- data.frame(B.s.glm)
SE.s.glm <- data.frame(SE.s.glm)

B.ns.glm <- data.frame(B.ns.glm)
SE.ns.glm <- data.frame(SE.ns.glm)


# Figure 2 Plot
windows(width=8, height=2.5)
par(mfrow=c(1,3))
n <- NROW(B.glm)

plot(0, 0, type='n', ylim=c(-.5, 1.75), xlim=c(min(anes$yr), max(anes$yr)),
  cex.lab=1.15, xaxt="n", yaxt="n", main="All Individuals", cex.main=1, xlab="Year", ylab="Income Coefficient")
  axis(side=1, at=c(1960, 1980, 2000), cex.axis=1.1)
  axis(side=2, at=c(0, 1), cex.axis=1.1 )
  abline(h=0, col="gray")
    for (i in 1:n) {
    points(p.yr[i], B.glm$z.inc[i], type='p', pch=19)
    segments(p.yr[i], (B.glm$z.inc[i]-SE.glm$z.inc[i]), p.yr[i],  (B.glm$z.inc[i]+SE.glm$z.inc[i]))
  }

plot(0, 0, type='n', ylim=c(-.5, 1.75), xlim=c(min(anes$yr), max(anes$yr)),
  cex.lab=1.15, xaxt="n", yaxt="n", main="Southerners", cex.main=1, xlab="Year", ylab="Income Coefficient")
  axis(side=1, at=c(1960, 1980, 2000), cex.axis=1.1)
  axis(side=2, at=c(0, 1), cex.axis=1.1 )
  abline(h=0, col="gray")
    for (i in 1:n) {
    points(p.yr[i], B.s.glm$z.inc[i], type='p', pch=19)
    segments(p.yr[i], (B.s.glm$z.inc[i]-SE.s.glm$z.inc[i]), p.yr[i],  (B.s.glm$z.inc[i]+SE.s.glm$z.inc[i]))
  }

plot(0, 0, type='n', ylim=c(-.5, 1.75), xlim=c(min(anes$yr), max(anes$yr)),
  cex.lab=1.15, xaxt="n", yaxt="n", main="Non-Southerners", cex.main=1, xlab="Year", ylab="Income Coefficient")
  axis(side=1, at=c(1960, 1980, 2000), cex.axis=1.1)
  axis(side=2, at=c(0, 1), cex.axis=1.1 )
  abline(h=0, col="gray")
    for (i in 1:n) {
    points(p.yr[i], B.ns.glm$z.inc[i], type='p', pch=19)
    segments(p.yr[i], (B.ns.glm$z.inc[i]-SE.ns.glm$z.inc[i]), p.yr[i],  (B.ns.glm$z.inc[i]+SE.ns.glm$z.inc[i]))
  }

# Figure 2
savePlot(filename="nes.yr", type=c("pdf"), device=dev.cur())
savePlot(filename="nes.yr", type=c("ps"), device=dev.cur())


